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Abstract. We present a parameter study of self-consistent models of protoplanetary disks around Herbig AeBe stars. We use 
the code developed by DuUemond and Dominik, which solves the 2D radiative transfer problem including an iteration for the 
vertical hydrostatic structure of the disk. This grid of models will be used for several studies on disk emission and mineralogy in 
followup papers. In this paper we take a first look on the new models, compare them with previous modeling attempts and focus 
on the effects of various parameters on the overall structure of the SED that leads to the classification of Herbig AeBe stars 
into two groups , with a flaring (group I) or self-shadowed (group II) SED. We find that the parameter of overriding importance 
to the SED is the total mass in grains smaller than 25yum, confirming the earlier results by Dullemond and Dominik. All other 
parameters studied have only minor influences, and will alter the SED type only in borderline cases. We find that there is no 
natural dichotomy between group I and II. From a modeling point of view, the transition is a continuous function of the small 
dust mass. We also show that moderate grain growth produces spectra with weak or no lO^m feature, both for flaring (Group 
I) and non-flaring (Group II) sources. The fact that sources with weak features have been found mostly in Group I sources is 
therefore surprising and must be due to observational biases or evolutionary effects. 
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> 1. Introduction 



Circumstellar Disks - circumstellar matter 



Herbig Ae/Be stars are young, intermediate-mass stars sur- 
rounde d by circumstellar di sks (see e.g. Waters & Waelkens 
d 1998b or iNatta et al.1 (I2OOOI) for a review). Solid-state mate- 
rial makes up about a percent of the disk mass, and can be 
studied in the thermal infrared. These studies indicate that the 
grains must have experienced both growth and thermal pro- 
cessing, since their spectra differ subst antially from that of in- 



terste llar dust ( Bouwman et al.. 200 ll: Ivan Boekel et al.ll2003l 



the buUc of the observable warm silicate grains in these sources 
have sizes larger than several microns. 

The general shape of the SED of Herbig stars has been a 
subject of study for a considerable time and is now believed 
to be quite well understood. The latest addition to this under- 
standing was the realization that the prominent emission close 
to 3fj.m, interpreted in earlier models interpreted as o riginating 
in accretion luminosity (e.g. lHiUenbrand et al ]ll992h . is due to 
a puffed-up inner rim of the disk. This rim is formed at a lo- 
cation given by the dust-evaporation temperature. Inside the 
rim, the disk does not contain any dust and is mostly trans- 
parent (optically thin) to the incident stellar light. The rim it- 
self c an absorb up to 25% of the stellar radiation d Natta et aT 



20051) . It is expected that this growth may lead ultimately to l200lh . This id ea was incorpora t ed into a more detailed phys- 



the formation of planets. The infrared spectra of Herbig stars 
can be divided into two groups (Meeus et al.,, ,2001). Though 
approximately similar in the near-IR, the far-IR flux of some 
is relatively strong (those of Group I) compared with that of 
others (Group II). This spectral difference is believed to reflect 
the overall shape of the disk. Group I sources have a geometri- 
cally flaring outer disk, which absorbs direct stellar light that is 
thenre -radiated at far-IR wavelengths. Group II sources have a 
flat geometry, which reprocesses far less stellar radiation in the 
outer disk. Many sources exhibit 10 micron silicate emission 
(these are labeled a, so either la or Ila). Some sources lack this 
emission (these are denoted lb or lib), possibly reflecting that 



ical model by Dull emond et al. I (1200 ih . who showed that this 
inner rim can indeed resolve the discrepancy between obser- 
vations and previous models in the wavelength region around 
3;um, with a self-consistent disk model. An important question 
is if this conclusion still holds in more complex models. 



Dullemond! (120021) and lDuUemond & DominikI (12004 here 



Send offprint requests to: J. Meijer, 



after DD04) developed far more sophisticated models includ- 
ing full 2D radiative transfer, coupled with self-consistent ver- 
tical structure calculations. These models confirmed that a self- 
consistent approach can explain many of the observed fea- 
tures of Herbig Ae SEDs. They demonstrated that the disk 
structure is closely linked to the properties and spatial dis- 
tribution of the dust grains. Although the vertical extent of 

the. disk is siinnnrtp.rl hv (las nrp.ssiirp. the. p.np.rsJV that hp.ats 
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the gas is obtained from the gr ains by dust-gas interactions 
(see iKamp & Dullemondl (120045 for details). The spatial dis- 
tribution of dust opacity that absorbs stellar light must there- 
fore control, to a large extent, the geometrical shape of proto- 
planetary disks. Since this opacity is provided mostly by grains 
smaller than about 25fim - that have large mass-extinction 
coefficients - this grain population is the mainly relevant one. 
Using this principle, DD04 performed a pilot study in which 
they showed that by increasing the total mass of the small 
grain component, one can alter the overall appearance of the 
spectral energy distribution from one qualitatively characteris- 
tic of Group II into Group I, and that this corresponds to either 
a non-flaring vs. a flaring geometry. A change in the geometry 
could also be achieved by a very steep power law in the surface 
density, which in eff'ect concentrates most of the disk mass in 
the inner rim, leaving the outer disk with very low mass and, 
consequently, opacity. 

This paper is the first in a series where we explore these new 
models in greater detail, building on the work of DD04. In this 
first paper, we focus on the general appearance of the disk 
spectra, and how it depends upon the various parameters en- 
tering the calculation. We present an extended parameter study 
of dust and disk properties of HAeBe star with the aim of bet- 
ter quantifying and understanding the cause of the Group I/II 
dichotomy more robustly as well as the presence or absence 
of the lOyum feature. We do not consider the detailed shape of 
the lOyum feature - mineralogical studies will be presented in 
further papers in this series. 

The model parameters that we address are - in addition to 
total disk mass and radial surface density gradient - disk size 
and inclination, some basic aspects of grain composition, and 
the grain-size distribution. To allow for a quantitative compar- 
ison between observations and models, we present the results 
in terms of a color-flux ratio comparison, first introduced by 
van Boekel et alJ (I20031 hereafter BWD03) to distinguish be- 
tween Group I and Group II sources. We then first study the 
SED of selected models directly, and later turn to a special 
color-color diagram in which we are able to display the results 
of the full parameter study. 

In Sect. |2l we describe the model assumptions, take a first 
look at the resulting SEDs, and present the basic model grid. 
In Sect. [3] we show the SEDs for a selected subset of our pa- 
rameter study. In Sect.|4] we describe the system to quantify the 
SED structure using a color-color diagram, and discuss in de- 
tail the parameter dependences of the model. We discuss some 
of the results in Sect. |5] and present our conclusions in Sect.|6l 



respect to the azimuthal angel <p vanish). The model consists of 
two parts: a 2D radiative transfer code (RADICAL/RADMC) 
and a ID vertical hydrostatic equilibrium code. The radiative 
transfer code calculates the dust temperature from the spatial 
distribution of the dust, and the dust opacity. The vertical hy- 
drostatic equilibrium code calculates the dust density from the 
dust temperature and vertical hydrostatic equilibrium. These 
two codes are applied alternately to an initial guess for the dust 
density distribution until convergence is reached. For details of 
the method we refer to the references given above. 

No scattering is currently considered by the models. Since at 
optical wavelengths the scattering opacity for silicate grains 
(thought to be a dominant dust species in proto-planetary 
disks) is approximately a factor of a few larger than the ab- 
sorption opacity (and is about equal to it at near-IR wave- 
lengths) one may expect that this may aff'e ct the disk geometry 
and em erging spectral energy distribution. lDullemond & Natta 
('2003') investigated this effect and concluded that because scat- 
tering reflects part of the stellar light away from the disk before 
it becomes thermalized, the disk becomes cooler at all vertical 
heights. However, the vertical surface height - defined to be the 
surface at which the stellar light is intercepted - will increase, 
simply because of the added (scattering) opacity component. 
The effect on the SED is that in the near- and mid-IR the flux 
may reduce by tens of percents, while at millimeter wavelength 
the reduction is a few percent at most. We conclude that ne- 
glecting the effects of scattering may cause a modest (system- 
atic) effect on the disk structure and spectral appearance in- 
feiTed, and therefore on the exact location of the Group I/II 
boundary. However, it does not affect selectively either Group 
I or II, and will therefore not affect any conclusions regarding 
the cause of the Group I/II dichotomy itself. 

2.2. Modifications for the present study 

The following modifications to the original code have been im- 
plemented for the computations presented in this paper: 



2.2.1 . Grain properties 

In the pilot study by DD04, the grain population was assumed 
to consist only of silicate particles with a radius a - 0. 1 fim. 
Here, w e account for both sih cate dLaor & Drainelll993b and 
carbon (IPreibisch et al.l Il993b grains that may range in size 



from a minimum value amin to a maximum value a^^x accord- 
ing to 



2. Model setup 

2.1. Tlie DD04 code 

To model the proto-planetary dis ks of Herb i g sys tems, 
we used the code develo ped by i Dullemond (2002) and 



DuUemond & DominikI (12004) in which 2-D radiative transfer 
is solved subject to the constraint of radiative equilibrium and 
vertical hydrostatic equilibrium for a passive disk (i.e. one in 
which the only source of energy is irradiation by the central 
star) for which the surface density structure is described by a 
power law. 

The code uses a spherical coordinate system (r, 0, </>) with the 



n(a) I 



(1) 



We ad opt m = 3.5, whic h is representative of interstellar 
grains (' Mathis et al. 1977 ). This power-law shape size dis- 
tribution is expected on theoretical gro unds whenever g rain- 
grain collisions may lead to shattering (IDohnanvi 1969 ). but 
differs from that expected on the basis of grain growth (e.g. 
Dullernond & Dominik 2005; Tanaka et al. 2005) in which 
case the dust size distribution has one or two peaks that 
shift with time. If the distribution results from an equilib- 
r ium between growth and fra gmentation, it may also be flat 
dDuUemond & DominikI 2005h . It is therefore not at all obvi- 
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choice. We used it nevertheless, because it is the standard as- 
sumption in most disk models. The relative abundance of car- 
bon and silicate grains is explored as a model parameter. 

Grains of different sizes are all assumed to have the same 
temperature. Technically this is achieved by averaging the 
mass absorption coefficients of the different grain sizes and 
materials, weighted by the relative mass contained in each 
grain size and material. This averaged mass absorption coef- 
ficient is then used to solve the radiative transfer problem. 

2.2.2. Small and large grains, and the small grain disk 
mass 

Many studies of dust evolution in disks discuss the presence of 
"large" grains. Depending on the type of study, this term may 
mean grains of a few fim, as they may be detected by studying 
the lOyum emission feature, or mm-sized grains as submm and 
mm studies report. It is therefore important to define clearly 
what is meant by "large" and "small" in the context of a study. 

In the present paper, we are interested mainly in the overall 
structure of the disk model, and this structure will be domi- 
nated by the grains that carry the opacity in the disk, i.e. the 
relatively small grains. Instead of using the total disk mass as 
an input parameter, we therefore use the term small grain mass 
of the disk as the one main parameter of our disk models. This 
covers all grains smaller than 25yum. In most models we only 
consider these grains and vary the small grain mass of the disk 
over a large range. This does not necessarily imply that the true 
dust mass in the disk equals to this value. In fact, one way to 
justify a large range in the small grain mass is by assuming 
that much of the original dust mass in the disk has been con- 
verted to larger grains, that no longer influence the structure 
of the disk. In Sect. 13.21 we demonstrate a computation of a 
few models in which we add significant amounts of dust mass 
in the form of a midplane layer of large grains, and prove that 
the general shape of the SED at near and mid-IR wavelength is 
largely unaffected by this mass. 

However, we study in detail the effects of changing the rel- 
ative abundance of grains below 25jum: these changes are ex- 
pected to affect both the strength of emission features close to 
lOjum, and the overall disk structure. 



2.3. First look: an example model and comparison to 
the DDN model 

To obtain an impression of the new self-consistent models and 
how they compare with the DDN models, we use the model fit 
completed bv Dominik et al.,(,2003i, hereafter DOB) for the case 
of AB Aurigae. We have not made a specific fit with the new 
model to the observed data. For the sake of comparison, we 
have taken the disk parameters exactly as they were derived by 
D03, and solved the radiative transfer problem and the vertical 
structure for this case. In Fig. [1] we show the observed SED 
of AB Aurigae and its D03 fit, the new model seen face-on, 
and the new model seen at an angle of 65 degrees (the value 
derived by D03). 

As can be seen, there are significant differences between 
the models, most important of which may be the difference 
in the 3//m bump, which represents the radiation that is caused 

nnr^otl\7 \w7 fhf^ inn <=»r" -nm W/liilfi \wt tVifi T~\0'5 fit tViio »-famr\n r\{- tVita 




Fig. 1. The observed SED of the Herbig star AB Aurigae 
with several model fits. The solid line shows the fit using the 
extended Chiang and Goldreich (IChiang & Goldreichl[l999l) 
- like model including the effects of a puffed-up inner rim 
(Dullemond et al. 2001). The fit was originally obtained by 
(iDominik et al.ll2003l) . using an inclination of 65°. The dotted 
line shows a face-on model obtained for the same star and disk 
parameters, the dashed line shows the same model tilted to a 
65°incIination. 



SED is well fitted, the new models underestimate the flux. This 
is a systematic effect: the self-consistent model s under-predict 



this fl ux systematically, a fact also noted by IVinkovic et al 



(2006). This is partly because the temperature drops rapidly 
into the rim, and part of the absorbed energy is emitted at wave- 
lengths far longer than 3 /urn. It may be in part because the ge- 
ometric shape of the inner rim is in reality more rounded than 
assumed here ( Isella & Nattal 20051) . This would also remove 
the extreme dependence of the strength of the inner rim emis- 
sion on the inclination of the disk, as can be seen by comparing 
the model spectra for a face on disk with one that is inclined 
by 65°. The unrealistically strong dependence on incUnation is 
also the reason why direct model fits infer high inclinations in 
out model - the inner rim emission is highest at these angles. 
Isella & Nattal ( 2005 ) demonstrated that a more rounded rim 
has its maximum flux in the face-on orientation. Unfortunately, 
we still need to devellop completely self-consistent models for 
the inner rim. We only need to notice that there is a problem, 
and that inclination fits should not be trusted. In Sect. I5.1l we 
will consider the systematic effects caused by the weak inner 
rim. 

Another side-effect of the model at 65°inclination is that the 
stellar flux is clearly reddened by the absorption of the outer- 
most disk parts. At this inclination, the flaring outer part of the 
disk extends to angles of up to 50°, so that the line-of-sight 
from the observer to the star passes through these outer parts. 
Dust grains present th ere absorb the stellar light and redden the 
star (see, for example I Whitney et al ] |2003h . 

The far-infrared wavelength range of the SEDs emits more 
flux when compared to the D03 model. This is a direct con- 
sequence of the different height of the inner rim. In the D03 
models, the higher rim shadows the outer regions of the disk 
more, and consequently less reprocessed radiation emerges 
from these regions. Clearly, it would be possible to construct a 

1 - 1 . i_ . : 1 „ 1 . . r. I A . 4.1 . : i i 



4 



J. Meijer et al.: A parameter study of self-consistent disk models around Herbig AeBe stars 



Table 1. The parameters describing the disk structure and grain 
properties. The power law describing the grain size distribution 
has a fixed index m = 3.5, representative of interstellar grains. 
The values of the parameters that are not varied are marked by 
a 



disk parameter 


adopted values 


log dust mass [Mq] 




-6.0 


-5.5 


-5.0 


-4.5 -4.0 


outer radius [AU] 




*200 


400 


600 




mass distribution 


P 


-1.0 


•-1.5 


-2.0 




dust parameter 






adopted values 


size range [pm] 




*0.01 


0.1 


1 






^max 


*1 


5 


10 


25 


mass fraction 


/sil 


•1 


0.95 


0.91 


0.83 




fcar 


•0 


0.05 


0.09 


0.17 



ing the small grain mass in the outer disk regions. We have not 
attempted this here. 

2.4. The model grid 

Our model grid was computed for a central star of mass 
M* - 2.5 Mq, radius - 2.1 Rq and eff'ective temperature 
Teff - 10 000 K, and therefore a luminosity = 39.5 L©. Its 
spectral energy distribution was represented by a Planck func- 
tion. An overview of the disk parameters is given in Table [T] 
The free parameters describing the global disk structure were: 
the total mass in small dust grains Mdust, the mass distribu- 
tion, and the outer disk radius /?disk- A fixed gas-to-dust ra- 
tio / = 100 relates M^ust to the total disk mass, although the 
amount of gas in the disk enters the calculation in no way. 
The radial distribution of the surface density was given by 
S(r) oc r^P, i.e. it is described by the power-law index p and 
the total disk mass. We note that although the outer disk ra- 
dius is a free parameter, the inner disk radius is computed self- 
consistently using the dust sublimation temperature Tsub- The 
dust is assumed to consist of silicate and carbon grains. In the 
current study, we were not interested in the effects of different 
sublimation temperatures of different species. We therefore as- 
sume that both species sublimate together at Tsub = 1 500 K, 
typical of silicates. We also investigated the effects of grain size 
and computed a set of models with different grain size ranges 
(see Table[T]i. The full model grid therefore comprises 50 mod- 
els. The relative abundance of the two dust species considered 
was also varied, but the effect was negligible. 

In Figs.|2]and|6l we show the main results of our parameter 
study. Both figures contain 8 panels each, and in each panel we 
study the dependence of the model spectra on one parameter. 



3. Spectral energy distributions 

Figure |2] shows the SEDs of the different models relative to our 
standard model, which is given by the solid curve in each panel. 
This model uses a total mass in small dust grains of IQ^^Mq, 
a surface density power law with p - -1.5, a dust size distri- 
bution between 0.01 and l.Qpm, a carbon/silicate mass ratio of 
and a luminosity of 4OL0. The SED of the standard model is 
computed at an inclination of 50°. In each panel, we then vary 



3.1. Overview over SED parameter dependence 

Figure |2] shows the spectral energy distributions for a subset 
of our parameter study. In each panel, we start from our stan- 
dard base model and then vary a single parameter The result- 
ing SEDs are shown, and the different values of the varied pa- 
rameter is shown in the panel legend. We present a far more 
detailed discussion of the parameter dependencies in Sect. 14. 31 
here, we only highlight the obvious properties of these calcu- 
lations. It is clear that the dependence of the SED on the mass 
of small dust grains (panel a), the inclination (panel b) and the 
maximum grain size (panel f) is far stronger than the influ- 
ence of all other parameters. In panels d, e, and h, hardly any 
variations are seen, whereas the data in panels c and g show 
moderate changes. 

One result is particularly significant: The shape and strength 
of the 3fj.m bump is extremely robust and similar in all mod- 
els, with only two exceptions. One exception is high inclina- 
tion models (/ > 75°). In these models, the inner rim is seen 
through the outer disk, and the 3fj.m flux is strongly absorbed. 
The other exception is the model with a luminosity of only 
IOLq for which the inner rim emission is less significant. The 
reason for this latter case is that at low luminosities, the dust 
evaporation radius is far closer to the star, and therefore deeper 
inside the gravitational potential of the star. The gravity of the 
star reduces the ratio of surface height to distance, so that a 
smaller fraction of the total luminosity is absorbed by the in- 
ner rim. Apart from those exceptions, this part of the spectrum 
shows an amazing lack of variation. This is a reflection of two 
important facts: (i) the distance of the rim is set by evaporation 
physics and therefore independent of the disk properties, (ii) 
its height is only weakly dependent on disk properties because 
the rim is strongly optically thick and remains so regardless of 
parameter changes. 

The SED variations in panel b (inclination) only occur be- 
cause of the location of the observer relative to the object. If 
we focus on the remaining parameters that alter the intrinsic 
properties of the star-disk system, the strongest influence is ex- 
erted by the mass in small grains, which can switch the outer 
disk from optically thick to optically thin. 

Some of the variations in panel f are due to a significant re- 
duction in the optical depth of the outer disk, which is caused 
by an increase in the upper limit of the dust grain size distribu- 
tion to 25/zm. In addition, the lO/im feature is also weakened 
because this feature is mostly due to grains with sizes below 
3/im. 

3.2. Large grains in the mid plane 



We considered a wide range of total small dust grain mass in 
the disk in order to produce models of significantly different 
properties. In reality, the amou nt of solids in a disk might spa n 
a far smaller range ov values (lAcke & van den Anckeiil2004l) . 
and the small dust mass required to reproduce group II models 
may be compensated for by larger grains or even pebbles and 
boulders in the disk midplane. DD04 already demonstrated that 
converting a large fraction of small grains into large midplane 
particles produces models that are similar to those of small 
dust mass, apart from the fact that the submm luminosity is 
increased. Therefore, throughout this paper, low dust mass can 
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X (fim) X (/if") 

Fig. 2. Panel representing an overview of the spectral energy distributions resulting from our parameter study. In each panel, 
we start from our standard base model (L* = 40Lq, Mdust = IO'^Mq, i = 50°, p = -1.5, Rdisk = 200AU, = O.Olyum, 
flmax = l//m) and vary a single parameter, as specified in the panels: (a) dust mass, (b) inclination, (c) index of power law of mass 
distribution, (d) disk size, (e) minimum grain size, (f) maximum grain size, (g) stellar luminosity, and (h) carbon/sihcate ratio. 



be assumed to correspond to models with an efficient conver- 
sion of micron-sized dust to far larger particles. 

To assess the effect of sedimentation on the spectral en- 
ergy distribution, we calculated an additional five models, in 
which we introduced a mid-plane layer of large grains. This 
was achieved by taking a model from the standard grid with 
a dust mass less than 5 x 10""^ M0 and adding a mid-plane 

Invpr nf 9 mm araine nf a maee eiirh that thp tntal HieV mnw 



equaled 5 x 10" M0. The resulting values of the mass in the 
midplane layer can be seen in Table |2] The large grain mass 
fraction varies between 0.998 for model MLO and 0.000 for 
model ML5, which is equivalent to its parent model M5 (see 
TableHi. 

The SEDs are shown in Fig.|3] In the upper panel, we com- 
pare all models that include a mid-plane layer. While there are 

etrnnff HifFprpnrpe in thp npar miH and far infrarpH thp mnHple 



6 



J. Meijer et al.: A parameter study of self-consistent disk models around Herbig AeBe stars 



Table 2. Dust mass in small (0.01 - 1 /im) and large (2 mm) 
grains for the mid-plane layer grid. The star has the following 
parameters: = 2 Mq, R, ^ 3Rq and T, = 10000 K, and 
thus L, = 80.7 Lq. The disk has a size of 200 AU and a dust 
mass of 5 X 10""^ M0. The small dust in the disk has a grain 
size range of 0.01 - 1 fim with a power-law distribution with a 
power of -3.5. The large dust in the mid-plane layer has a size 
of 2 mm. 



-Mdust.orig 

model [Mq] 



MLO 
MLl 
ML2 
ML3 
ML4 
MLS 



1 X IQ-^ 
3 X 10-'' 
1 X 10^^ 
3 X 10-5 
1 X 10-* 
5 X 10-^ 



^'^midplane 

[Me] 
4.99 X lO-* 
4.97 X 10-" 
4.90 X 10-" 
4.70 X 10-" 
4.00 X 10-" 




large grain 
mass fraction 



0.998 
0.994 
0.980 
0.940 
0.800 
0.000 















ML5 = M5 

ML4 

- - - ML3 

: WL2 

ML1 

MLO 





10.0 
X (/urn) 




Fig. 3. The SEDs of the models in which a midplane layer of 
large grains is added so that the total dust mass in the disk is 
kept constant. The upper panel shows the models MLO. . . MLS 
that display a transition from group I to group II SED type, but 
all converge again at about 800 fim. The lower panel compare 
each pair of models with equal small grain mass. The SED's 
are shifted vertically for improved visibility. 



converge in the submm region, indicating equal total mass. In 



small grains, one without the added midplane layer, one with 
the added layer. It is clear from these plots that the shape of 
the SED blueward of about 60 yum is determinated by the mass 
in small grains, while at wavelengths longer than 100 /vm clear 
differences can be seen. In a classification of SEDs based on 
the SED up to 60//m, we can therefore ignore the presence 
or absence of a mass-compensating midplane layer of large 
grains. We also demonstrate this in Sect. 14.21 

This implies that sources with the same small grain mass 
have the same group classification, but not necessarily the 
same total dust mass because mass may be 'hidden' in 
large grains. Since large grains affect the shape only at 
sub-mm wavelengths a nd are therefore not evident in the 
van Boekel et al.l (l2003h diagram, adding large grains in the 
mid plane of the disk does not influence the Meeus classifica- 
tion significantly. We can therefore conclude that differences 
between the properties fo group I and a group II sources are 
due to the small (a < 1 fj.m) grain mass. While the models in- 
fer a factor of 100 difference in sub-mm flux between group 
I and II, observations only show a factor of 10. This implies 
that large grains must exist in the group II sources to increase 
the sub-mm flux to the observed levels. As demonstrated by 
Acke et al.l (120041) . the sub-mm slope of group II sources also 
indicates that there must be more large grains present than in 
group I sources. 

4. Color-color diagram 

4.1. Quantifying SED classification 

To quantify the amount of flaring in all calculated models, 
we needed to employ a method that is both fast and accu- 
rate. Previously, the classification of group I/II was achieved 
by fitting power-la w and/or black-body functions to the SED 
(i Meeus et al.ll20011) . splitting the SEDs into two groups based 
on the presence (group I) or absence (group II) of a black-body 
component in addition to the power-law component. Since this 
proved to be an inefficient method when analyzing hundreds of 
SEDs, we chose a simple quantitative way of classifying SEDs. 
BWD03 provided us with such a method. 

In BWD03, an IRAS - OT6o color versus Lnir/Lfir dia- 
gram was presented. Lnir was the integrated luminosity from 
JHKLM photometry and LpjR the same quantity derived from 
the IRAS 12, 25, and 60 //m data points. Figure |4] shows these 
two luminosities for two sources typical of group I and group 
II. The IRAS mi2—m^ color is sensitive to the slope of the SED 
between 12 and 60 yum, which is steeper in group II sources re- 
sulting in a lower IRAS nin - wJeo color for group II sources 
than group I sources. The near to far infrared luminosity ratio 
Lnir/^fir is sensitive to the disk surface-temperature distri- 
bution. Since group I sources have warmer outer disk surface 
temperatures than group II sources their FIR emission will be 
higher relative to the NIR emission and the NIR over FIR ratio 
will be lower. Data for group II sources occupy mainly the up- 
per left corner whereas group I sources will occupy primarily 
the lower right corner An example of the diagram is shown in 
Fig.|5] The sources used in Fig.|5]are listed in Table[3] 

The line Lnir/^^fir - (mn - 'weo) + 0.9 was determined 
to be the most reliable means of separating group I from 
group II sources. Sources above the line belong to group II, 
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2: HD 100546 (lo) 




Fig. 4. SEDs of two sources HD 100546 and HD 163296 show- 
ing the integrated fluxes Lnir and Lfir- The ratio of these inte- 
grated fluxes and the IRAS - color is used in all van 
Boekel diagrams in this paper. 
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Fig. 5. The Ivan Boekel et"al (2003) diagram for the sources 
listed in Table[3] On the ordinate, the IRAS m\2-ma) color and, 
on the abscissa, the LNiR/i^FiR fux ratio is given. The num- 
bers of the observations correspond to the numbers in Table [3j 
Pluses indicate la sources, crosses lb sources, diamonds Ila 
sources, and the square a group lib source. We note that all lb 
sources are in the upper right corner of the group 1 zone. The 
dashed line represents LNiR/ipiR = (pin - meo) + 0.9. 



doe s not comply w ith this r ule because of its di fferent geome- 
try ( Fukagawa et al . 2006; Fuiiwara et alj|200"6l) . However, the 
porperties of inner part of its disk indicates that it is group 11 
jLeinert et aDl2004 . 

Other aspects of the disk and the dust also play a role in 
determining the position of a model in this diagram. A model 
without a puffed-up inner rim emits less NIR and more FIR 
flux than a model with a puffed-up inner rim. This is because 
the puffed up inner rim re-emits more radiation in the NIR and 



FIR flux will decrease more on the blue side and increase in 
slope because the shadow does not affect the outer parts of the 
disk in which the flux redward of 20 /vm originates. Therefore, 
a model without a puffed-up inner rim has a lower NIR/FIR 
ratio and a redder IRAS m\2 - mgo color and will be shifted 
down and to the right in the diagram with respect to the models 
with a puffed-up inner rim. 

Adding a mid-plane layer of large grains to a model will 
increase the FIR flux longward of approximately 50 //m. The 
large grains in the mid plane are only visible in the outer parts 
of the disk, where they are cold and emit at longer wavelengths. 
This means that adding a mid-plane layer will increase both the 
FIR flux and the IRAS nin - meo color, shifting the predicted 
data points for a model with an added mid-plane layer of large 
grains to the right of and slightly lower than the original model. 

The NIR and FIR windows were chosen so that there is only 
a limited contribution from the sihcate emission around 10 /im 
and it is primarily the geometry of the disk determine the po- 
sition in the diagram. However, (cold) crystalline silicates also 
emit a broader and less conspicuous feature that peaks around 
20 jim. If there is a large amount of crystalline silicates in 
the outer pa rt of the disk (beyond 2 AU), such as that in HD 
142527 (see Ivan Boekel et al.ll2004l) . this will contribute sig- 
nificantly to the FIR flux. This means that sources with highly 
crystalline (outer) disks will appear lower in the diagram. 

Figure |6] shows the position of all models in our parameter 
study in the diagnostic color-color diagram. We varied each 
parameter for 6 different small-grain dust masses, using the 
models M0-M5 as the basis for each variation track. 



4.2. Influence of large mldplane layer grains 

In Fig. |2] we show the positions in the diagnostic diagram of 
the models M0-M5, and compare them with the positions of 
the corresponding models in which a low value of the small 
grain mass is compensated by a layer of large grains in the 
mid-plane. This figure clearly confirms our expectation from 
Sect 13.21 The positions of the models with the same small- 
grain mass are very close to each other. Models ML2-ML5 
occupy exactly the same position as the corresponding models 
M2-M5. A visible shift occurs only for models with the low- 
est small grain masses, MLO/MO and MLl/Ml. However, the 
shift is small compared with the overall range for both observa- 
tions and modelhng results in this diagram. We may therefore 
restrict our discussion to the small grain masses. 

4.3. General parameter trends 

We will discuss the trends observable in each panel of Fig |6] 
and the physics in the models responsible for these trends. 

Dust mass In Figs. |2ja) and|6ja), the effect on the spectrum 
of increasing the dust mass can be seen. As concluded be- 
fore in DD04, dust mass is an important factor in determin- 
ing the disk geometry. Higher mass provides hgher opac- 
ity, which causes more flaring. An increasing dust mass 
therefore shifts a model from group 11 to group 1. It should 
be noted, however, that primarily small grains (<10 fj.m) 
contribute to the opacity in the disk, so the dust mass men- 
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Table 3. A list of sources that were used to compare with our models. In Col. 1 we give the identification number that appears in 
all diagrams. Column 2 contains the name of the star . In Col. 3, 4, 5, and 6, th e distance, effective temperature, stellar luminosity 
and stellar mass are given. This data was taken from lvan Boekel et al.l (l2005b . In Col. 7, the stellar radius is given. In Col. 8, we 
indicate the IRAS mi2 - »i6o color. Column 9 lists the flux-ratio of NIR over FIR. The NIR flux is determined from JHKLM 
photometry. The FIR flux is derived from IRAS 12, 25 and 6 fluxes. In Col. 10 a nd 1 1 the ratios of NIR and FIR flux to stellar 
flux are given. In Col. 12 the group classification according to lMeeus et al.l (1200 ih is given. 
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(4) 


(5) 


(6) 
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star 


d 
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M 


R 


[12] - [60] 




^nir/^* 




group 






[pc] 


[K] 


[Lq] 


[M©] 


[R©] 












1 


AB Aur 


144 


9528 


47.0 


2.4 


2.5 


1.504 


1.763 


0.274 


0.155 


la 


2 


HD 100546 


103 


10495 


32.0 


24 


1.7 


1.000 


0.281 


0.113 


0.404 


la 


3 


HD 179218 


243 


10495 


100.0 


2.9 


3.0 


0.460 


0.638 


0.099 


0.155 


la 


4 


HD 97048 


175 


10000 


43.7 


2.5 


2.2 


1.806 


0.587 


0.095 


0.162 


lb 


5 


HD 100453 


111 


7396 


7.9 


1.7 


1.7 


1.840 


1.089 


0.287 


0.263 


lb 


6 


HD 135344 


140 


6592 


8.1 


1.6 


2.2 


2.943 


3.545 


0.360 


0.102 


lb 


7 


HD 139614 


140 


7852 


8.1 


1.7 


1.5 


1.649 


0.897 


0.191 


0.213 


lb 


8 


HD 169142 


145 


8204 


14.5 


2.0 


1.9 


2.503 


1.026 


0.136 


0.133 


lb 


9 


HD 104237 


116 


8414 


35.0 


2.3 


2.8 


-0.492 


3.599 


0.252 


0.070 


Ila 


10 


HD 142527 


198 


6252 


29.0 


2.5 


4.6 


2.513 


1.844 


0.417 


0.226 


Ila 


11 


HD 142666 


145 


7580 


13.5 


1.8 


2.1 


-0.162 


1.761 


0.222 


0.126 


Ila 


12 


HD 144432 


145 


7345 


10.2 


1.8 


2.0 


-0.291 


2.319 


0.324 


0.140 


lla 


13 


HD 150193 


150 


8974 


24.0 


2.3 


2.0 


-0.957 


1.654 


0.212 


0.128 


lla 


14 


HD 163296 


122 


8730 


24.0 


2.0 


2.2 


0.763 


2.827 


0.286 


0.101 


lla 


15 


51 Oph 


131 


10000 


245.0 


2.0 


5.1 


-2.924 


9.294 


0.066 


0.007 


lla 


16 


HD 95881 


118 


8990 


11.6 


1.7 


1.4 


-1.971 


2.719 


0.203 


0.074 


lib 



Table 4. Flux ratios and infrared colors for selected models. The first sections contains the different AB Aur models shown in 
FiglH The second section shows the values for the standard models with diff'erent small grain masses. The final section shows 
values derived from the models with a midplane layer of large grains. 



model Mdust i [12]-[60] Fnir/Ffir Fnir/F. Ffir/F. 

Me 



CGplus 


5(- 


-2) 


65 


2.158 


1.680 


0.273 


0.162 


RADICAL 


K- 


-3) 


6 


1.793 


0.401 


0.142 


0.355 


RADICAL 


K- 


-3) 


65 


1.671 


0.802 


0.206 


0.256 


MO 


K- 


-6) 


50 


0.08 


1.60 


0.156 


0.097 


Ml 


3(- 


-6) 


50 


0.61 


1.26 


0.164 


0.130 


M2 


K- 


-5) 


50 


1.01 


1.00 


0.167 


0.167 


M3 


3(- 


-5) 


50 


1.24 


0.86 


0.173 


0.202 


M4 


K- 


-4) 


50 


1.41 


0.74 


0.177 


0.239 


M5 


5(- 


-4) 


50 


1.34 


0.63 


0.190 


0.303 


MLO 


5(- 


-4) 


50 


0.42 


1.69 


0.154 


0.091 


MLl 


5(- 


-4) 


50 


0.72 


1.27 


0.161 


0.127 


ML2 


5(- 


-4) 


50 


1.01 


0.99 


0.165 


0.166 


ML3 


5(- 


-4) 


50 


1.31 


0.89 


0.173 


0.194 


ML4 


5(- 


-4) 


50 


1.42 


0.74 


0.177 


0.238 



when we add a mid-plane layer of large (2 mm) grains in 
Fig. 13 

Inclination In Figs.|2|b) and|6jb), the effect of changing the 
inclination on the spectrum of models with different dust 
mass is shown. Increasing the inclination from face-on to 
the angle at which Ay equals 0. 1 (indicated by a triangle) 
enables more of the vertical inner rim to become visible. 
This causes the NIR flux to increase. Furthermore, increas- 
ing the inclination caused the outer part of the disk to ab- 
sorb an increasing amount of the inner rim radiation. This 
decreases the NIR flux, and reddens the entire spectrum 
making the IRAS mi2 - ma) color increase. At extreme in- 
clination angles (> 80°, not shown here), the reddening 
becomes so severe that the SED would no longer be classi- 



fied as a (Lada-)class (lLadalll987h II source, but rather as a 
class I source. 

Power-law of [Mass Distribution In Figs. |3c) and|6tc), we 
can observe the effect of changing the surface density 
power law. Steepening the power law while keeping the 
disk mass the same, causes relatively more of the mass 
to be on the inside of the disk. Between the shallowest 
{p - -\ .0) and the steepest {p - -2.0) power law slope, 
the mass in the inner 1 AU increases by a factor of 46 and 
the height of the inner rim by a factor of two. Having rel- 
atively more mass on the inside implies that there is less 
mass on the outside where the FIR flux originates. Thus 
the FIR flux decreases while the NIR flux increases and 
the NIR/FIR ratio increases as the power-law slope steep- 
ens. In the higher mass models the mass in the outer disk 
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Fig. 6. Plot with Ivan Boekel et al. I (l2003h diagrams, for variations of the most important parameters. In each panel, the dashed 
line represents LmR/LpiR = ('«12 - ^eo) + 0.9, the line used empirically to classify sources. The standard models which are shown 
in the first panel have an inclination of 50°, a surface density distribution power law of p = -1 .5, an outer radius of 200 AU and 
a grain-size distribution ranging from 0.01 to 1 fim. The parameter varied in each panel is: a) dust mass Dust mass increases to 
the right. The grey line connects the five standard models and one optically thin model that lies outside the diagram. This line 
is present in all panels, b) Inclination The standard models are at the average inclination of 50°. Inclination increases along the 
curve from the bottom to the right. The triangle in each curve indicates an Ay of 0. 1 mag, the point when the outer disk starts to 
significantly absorb the radiation coming from the inner rim. c) index of power law of mass distribution The power law becomes 
shallower from top to bottom, d) disk size, e) minimum grain size.f) maximum grain size, g) Stellar luminosity, h) Carbon-Silicate 
ratio. 



region is colder and emits mainly redward of 20 pm. The 
decrease in FIR flux blue-shifts the IRAS mn - meo color 
value. 

Outer Radius In Figs.|2|d) and|6td), the effect of an increas- 
ing outer disk radius is shown. Increasing the disk size 
while keeping the dust mass the same causes more of the 
mass to he on the outside and less on the inside. However. 



the mass is redistributed far less significantly than when 
the power-law slow is altered. The mass in the inner AU 
reduces by only a factor of two and the mass on the outside 
shows no significant increase. To achieve a similar redistri- 
bution of mass as in the power-law models, we would have 
to increase the disk size to 4 x 10^ AU. The position of the 
models in the. Hiauram would then hardlv chanue 
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Fig. 7. IRAS mi2 - ma) colour versus Lnir/^fir diagram with 
two different grids of models. Models MLO to MLS contain a 
mid-plane layer of large grains. The properties of this layer can 
be seen in table |2] Models M0-M5 are the original models, to 
which the mid-plane layer was added. 



Fig. 8. Effect of the strength of the inner rim on the position 
of a source in the color-color diagram. The diagram shows the 
locations of the measured SED of AB Aurigae, of the D03 fit 
using the DDN model, and of face-on and 65° SEDs of the new 
model, computed for the same disk parameters. 



Minimum grain size In Figs. |2e) andl^e), the effect of in- 
creasing the lower limit on of the grain-size distribution is 
shown. Because of the m - -3.5 power-law slope in the 
size distribution we use, most mass is in the largest grains. 
Increasing the minimum grain size therefore has little effect 
on the SED. The small grains emit mostly in the 10-20 //m 
wavelength range. This affects the FIR flux and the 12 jim 
IRAS color, and the NIR/FIR and the IRAS -meo color 
therefore increase. 

Maximum grain size In Figs. 13 f) and|3f), the effect of in- 
creasing the upper limit on the grain-size distribution is 
shown. By increasing the upper boundary of the size 
distribution power law while keeping the lower bound- 
ary constant, mass is extracted from the population of 
smaller grains and placed into the larger ones. Large grains 
have lower opacity (absorption cross section per mass). 
Consequently, the overall IR flux decreases and the lOyum 
and 20//m features reduce in strength. The FIR flux de- 
creases more rapidly than the NIR flux because the outer 
parts of the disk become optically thin (leading to weaker 
flaring), while the inner parts hardly change in structure. 
Beyond 60 fim, the flux increases as more and more cold 
material is revealed. However, when the maximum grain 
size exceeds 10 ;um, the flux also decreases in this part of 
the SED because the decrease in opacity becomes signifi- 
cant. This effect causes the reversal in the IRAS m\2 - m^t 
color. 

Stellar luminosity In Figs. |2jg) and|6fg), we show the effect 
of varying the stellar luminosity. Changing the stellar lu- 
minosity influences the spectra in two ways. The inner rim 
becomes much less pronounced in the calculations for a 
star with - lOL©, an effect already noted by DDNOl. 
At the same time, the low-luminosity object contains more 
flux at long wavelengths, between 60 and 100//m. In partic- 
ular, this effect leads to a horizontal shift in the diagnostic 
color-color plot. 

Carbon/silicate ratio In Figs. |2fh) and Hth), the effect of 
varying the carbon/silicon mass ratio in the dust opacity 
calculations is shown. It is immediately clear that this ef- 

j: i : ] 1 „£c A„ 4.1 1 „i _r 



SED. We will discuss in a companion paper (Meijer et al 
2007, submitted to A&A) how this is an important param- 
eter in terms of the strength of the lOjum feature. However, 
for the SED shape as probed by our diagnostic diagram, 
the carbon content has little or no effect. Since the addition 
of carbon increases the continuum level, both the NIR and 
FIR flux increase at the same rate, so that the ratio does not 
change. Since the slope between 12 fim and 60 fim hardly 
changes, also the IRAS m\2-m(,Q color remains unchanged. 

5. Discussion 

5.1. Effects of the weak inner rim 

In Sect. 123] we demonstrated that the inner rim emission in the 
new models is systematically weaker than in the simpler DDN 
models. In TableH] we quantify this result. Looking at the flux 
ratios derived from observed SEDs (see table |3]l, we can see 
that the rim emission (the near-IR flux, listed in Col. 10) ac- 
counts for the reprocessing of typically 20-30% of the stel- 
lar luminosity. In the models, the typical value is 15-17% and 
varies only weakly with disk properties. Clearly, this should 
also have a systematic effect on the diagnostic diagram. In 
figlH the infrared colors of AB Aur are compared with model 
predictions. Clearly, the difference is mainly a shift in the ratio 
of NIR to FIR luminosity in the disk. Taking this into account, 
it is unsurprising that our model predictions are below the data 
points of the observations. 

5.2. Effects of moderate dust growtli 

An important part of our parameter study has been considering 
the effects of moderate dust growth, i.e. either extending the 
upper limit of the size distribution from Ifim to values of up 
to 25/im, or the complete removal of small grains, modeled 
by shifting up the lower boundary of the size distribution. In 
Sect. 14.31 we found that removing grains below l/um has only 
a small effect on the structure of the models, while shifting 
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Fig. 9. All models in the entire parameter study that for- 
mally show a weak lOyum feature. We measured the peak- 
to-continuum ratio by drawing a linear continuum between 8 
and 13 jum and measuring t he peak flux in the way described 
by Ivan Boekel et alJ (l2003h . Pluses, crosses, diamonds, and 
squares show the location of observed sources of group la, lb, 
11a, and lib, respectively. Stars indicate models with all grains 
sizes below 1 yum removed. Triangles are extreme inclinations 
with weak features caused by self-absorption. 



removing smaller grains as well, including grains with sizes 
close to l//m) may affect the global SED. In this section, we 
study the basic effects of changes in grain size on the strength 
of the lOjum feature. Again, we do this without paying attention 
to mineralogical information - only the overall strength of the 
feature is considered. 



5.2.1. Weak 10//m feature 

Dust emission features in disks are an excellent way of prob- 
ing the composition and size of grains present in the disk. The 
observed IQ^ m features can vary strongly, both in strength 
and shape (e.g. Ivan Boekel et^al.ll2003l:IPrzygoddaet alJ 



van Boekel et alJbOOSHKessler-Silacci et alJl2006l) . In particu 
lar, there are a few sources that show no di scernible 10 yjm fea- 
ture at all. These objects were classified bv lMeeus et al. I (l200l]) 
into sub class b. O bservationally, at least in the bright ISO sam- 
ple (Me eus et alJ [2001). a significant fraction of objects with 
SED group lb were found, but no clear cases of group lib. 
An important question is, whether this is an automatic conse- 
quence of the disk flaring geometry (i.e. group 1 versus group 
11), or if additional factors are needed to produce weak features. 
We therefore measured the feature strength in all our models. 
Most of the models show a strong 10//m featu re and would 



therefore be classified as group la or 11a in the iMeeus et al 



scheme. However, there are several corners of parame- 
ter space where a weak silicate feature can exist. These sources 
would then have to be classified as group lb or lib. An overview 
plot of all models with a feature-t-continuum ratio lower than 
1 .05 (see Fig.|9]l. There are two basic ways of producing a weak 
emission feature: 

1 . Viewing a disk model at an inclination close to but not fully 
edge-on. In this case, the emission feature produced by the 
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Fig. 10. Example SEDs of two models with weak or absent 
10 fim feature. The group 1 model has a total dust mass of 
3 X 10"^ Mo, the grouop II model has 3 x 10 Mq. The other 
parameters are the same for both models: Rout = 200 AU, 20% 
carbon content, p - -1.5, and a size distribution for the dust 
grain of between 1 and lOyum with a slope m - -3.5. 



further out in the disk. A typical inclination for this eff'ect 
to occur is 65°. It is clear that this requires fine-tuning of 
the inclination in the modeling. Furthermore, models with 
a self-absorbed feature have significantly higher redden- 
ing. Observed sources with classification lb do not exhibit 
higher reddening than group la sources. Therefore, we are 
left with the second possibility: 
2. Efficient removal of grains smaller than about S/um . In the 
framework of our model, this can be achieved by either in- 
creasing the lower limit of the grain size distribution to a 
few yum, or increasing the upper limit of the grain size dis- 
tribution. If the grain size distribution is sufficiently steep 
(i.e. m < -3.5), most of the mass resides in large grains, so 
that increasing the upper limit effectively reduces the mass 
in small grains. We show two example of SEDs without 
strong features in FigfTOl 

While possibility (2) providesdelivers a reasonable expla- 
nation for the existence of the subgroup b sources, an im- 
portant problem is introduced by the fact that subgroup b 
appears to be much more abundant in group I sourc es than 
in group II sources (' Acke & van den Anckei 2004), while 
the models show just as many sources that could be clas- 
sified as lib. Therefore, we must be dealing here with 
eithe r an obs ervational bias, or with an evolutionary ef- 
fect. lAcke & va n den Ancker (2004) suggested that group Ila 
sources are systematically brighter than group lib sources, so 
that observational biases could be at work. 



6. Conclusions 

We arrive at the following conclusions. 

1 . The Parameter range investigated covers most of ob- 
served SEDs 

By comparing all models with observations, we are able to 
explain most, but we are unable to account for teh observa- 
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(11), HD 135344 (6) and HD 169142 (8). These stars all 
show evidence of disk gaps in imaging observations, which 
may be part of the reason why they are more difficult ob- 
jects to describe. 

2. The mass of small grains largely determines the SED 
type 

We find that the mass in small grains is the single most 
important parameter determining the SED type. While ex- 
treme values of the surface density power law (DD04) can 
also shift a disk from group 1 to 11 and vice versa, the val- 
ues required for this to occur seem too extreme. All other 
parameters have very limited influence on the classification 
of a source as group 1 or group 11. 

3. No group I/II dichotomy 

From a modeling point of view, we find that there is no 
clear dichotomy between group 1 and group 11 models. The 
position of a model in the diagnostic color-color diagram 
changes smoothly across the dividing line as parameters 
vary. The line connecting models of diff'erent small dust 
mass is almost perpendicular to the dividing line, while 
changes in other parameters cause shifts that are at least 
partially parallel to the dividing line. As the dust mass is 
changed in logarithmic steps, there is a gradual transition 
from fully flared (group 1) to non-flaring (group 11) disks. 
Some observed sources are in fact close to the dividing line, 
in particular AB Aur (source number 1) and HD 163296 ( 
source number 14) fall into this class. A scarcity of inter- 
mediate objects, if confirmed in studies of larger samples, 
would indicate a rapid transition if the distinction is inter- 
preted as an evolutionary effect. 

4. The strength of the inner rim emission is hardly depen- 
dent on disk mass 

The emission caused by the inner rim is hardly dependent 
on model parameters, as long as the dominant grain size is 
the inner rim remains the same. As long as the inner rim is 
sufficiently optically thick, its height depends only weakly 
on the details. In particular, the small grain mass influences 
only the inner rim emission when the surface density is de- 
creased by more than a factor of 100, compared with the 
case in which all dust mass is in small grains. This is con- 
sistent with observations tiiat show the 3A<m bump to be 
similar in all Herbig stars ( Meeus et alJl2001 ). 

5. The 3yum bump in the predicted spectra is too weak 
The models show that with the current assumptions, the 
self-consistent models do still not produce a sufficient 
amount of emission the 3//m bump. This causes a sys- 
tematic shift of the models in the color-color diagram. A 
more detailed implementation of the inner rim in models is 
highly desired. 
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